Time evolution in the Morse potential using supersymmetry: 
dissociation of the NO molecule 
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We present an algebraic method for treating molecular vibrations in the Morse potential per- 
turbed by an external laser field. By the help of a complete and normalizable basis we transform 
the Schrodinger equation into a system of coupled ordinary differential equations. We apply our 
method to calculate the dissociation probability of the NO molecule excited by chirped laser pulses. 
The dependence of the molecular dipole-moment on the interatomic separation is determined by a 
quantum-chemical method, and the corresponding transition dipole moments are given by approx- 
imate analytic expressions. These turn out to be very small between neighboring stationary states 
around the vibrational quantum number m = 42, therefore we propose to use additional pulses in 
order to skip this trapping state, and to obtain a reasonable dissociation probability. 
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An important model for anharmonic molecular vibra- 
tions, allowing also for dissociation is the Morse Hamil- 
tonian [uj 
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We shall consider the time development of a molecular 
state in this potential coupled also to an external classical 
field via the dipole interaction. Thus we shall give the 
solution of the time dependent Schrodinger equation: 



ih-^ = [H -$E(t)}\$), 



(2) 



where fi(X) is the molecular dipole moment, and E is the 
electric field strength of the external field. Introducing 



and the dimensionless operators X = aX, P = P/ah 
the dimensionless Hamiltonian Ho = {(2s + l)/hu>o}Ho 
can be recast in the form 



H = A*(s)A(s) - s 2 



(3) 



Here A(q) = ql - (s + 1/2) exp(-X) + iP, and A\q) = 
ql — (s + 1/2) exp(— X) — iP are generalized supersym- 
metric ladder operators 0], allowing the case q ^ s, as 
well. They satisfy the following commutation relations 

[A{q),A\q')\ = (q + q')I-(A(q) + A\ q >)), (4) 
[A(q),A(q')} = 0, [At(«),AV)] =0. 

As A^(s)A(s) — A(s)A^(s) — 2s+l is a shape invariant su- 
persymmetric operator Q , the solution of the eigenvalue 
equation Ho \ip m {s)) — E m (s) \tp m (s)} can be found as 
follows [g 



where AC, 



(2s-2r. 



I(2s-m)! ' and l^«( S ~ m )) is the sin S le 

nondegenerate vector annihilated by A(s — m) 



A(s — m) \ipo(s — m)) = 0. 



(G) 



Using coordinate representation it can be simply proven 
that |V>o( s — m ))i as we U as iV'm('S)} are normalizable if 
and only if m < s, therefore the number of bound states 
is [s] + 1, where [s] denotes the largest integer that is 
smaller than s. The infinite dimensional Hilbert space 
of the problem is, however a direct sum: H = HT ©7i + , 
where TC~ is the finite ([s] + 1) dimensional subspace 
spanned by the bound states, while H. + is an infinite 
dimensional subspace with elements obtainable as con- 
tinuous superpositions of positive energy eigenstates. 

Instead of using the non-normalizable continuum en- 
ergy eigenstates, we introduce algebraically a true or- 
thonormal basis in Ji that allows a natural discretization 
of time dependent problems. Let a = s — [s], and starting 
from the state \<fio) defined by the relation A(a) \cf>o) = 0, 
we introduce the following series of states 



\tn) = {l[C- 1 Al(a + k-l)}\<j>o), 



(7) 



k=l 



where the coefficients Ck = \Jk(k + 2a — 1) ensure the 
normalization. We will call the states \4> n ) quasi- number 
states Q . From the definition (Q) it is also seen that the 
neighboring quasi-number states are connected as 



A^(a + n) \4> n ) = C n+ i\4> n+ i) 
A(a + n)\^ n ) = C n \cj>n-i). 



The orthogonality of the \<p n ) states can be seen by using 
the commutators (||). In coordinate representation they 
have the form 



hM*)) = J^ m A\s)---A^s-m + l)\Ms-m)), = V^/T(2a + n)y a {x) exp( — ) 



y{ X )^ L 2a-i 
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-(s — to) 2 , to = 0, 1, ... [s], 



(5) 



(vW), 

(9) 



2 



where L 2 ^ -1 is a generalized Laguerre polynomial || |(| 
of the variable y(x) = (2s + l)e~ x . In contrast to the 
energy eigenfunctions, these wave functions form a com- 
plete orthonormal set of square integrable functions 
therefore the quasi-number states constitute a true or- 
thonormal basis in the Hilbert space. 

Using the properties of the supersymmetric ladder op- 
erators and A, one can calculate the matrix of Hq in 
this basis 

<&»| H \<j> n ) = (C m - s 2 + (m - [s]) 2 ) S min 
+ {[s] - n)C m 6 

m,n+l 

+ ([a] - m)C n S 

m+l,n ■ 

(10) 

It consists of two blocks, the one corresponding to 
< 7i, m < [s] is [s] + 1 dimensional, whereas the other 
one, corresponding to n, m > [s] is infinite dimensional, 
and both blocks are tridiagonal. The block structure 
shows that the first [s] + 1 quasi-number states span an 
invariant subspace of Hq, in which it can be diagonal- 
ized. This subspace should coincide with H.~ , as Hq has 
exactly [s] + 1 nondegenerate eigenvalues with normaliz- 
able eigenstates. The orthogonality relation 

(0n| il>m) =0, for m = 0, 1 . . . [s] , and n > [s] (11) 

equivalent to the above statement can be proven directly, 
and the expressions of the overlaps (</>„ ip m ) for n < [s] 
can be determined explicily . 

We will consider the time evolution of molecular vibra- 
tions subject to an external classical field as described by 
Eq. (0). In the dimensionless units used in Eq. (|J) the 
interaction term, H{ nt = —p(X)E(t) can be written as 



Hi, 



(12) 



where we have introduced ju, = ^-fl, with an effective 



charge q e = ^l\ x=0 , and E(t) = &±£gs.E(t) denotes 



the electric field strength in our units. As we will solve 
the Schrodinger equation (0) in the quasi-number state 
basis, we also need to know the matrix of the dipole mo- 
ment in this basis. Due to the algebraic properties of 
the quasi-number states, the matrix of n(X) can also be 
determined analytically for certain cases. Let us suppose 
that the dipole moment can be approximated by an op- 
erator of the form 



(13) 



where <Zj, di and ji are real numbers. Taking only a 
single term with d = 7 = would mean the simplest 
linear dipole moment borrowed from atomic calculations 
and used sometimes for molecules as well. This latter 
approach, however, overestimates the strength of the in- 
teraction at larger atomic separations and leads to an 
unrealistically high dissociation probability. 

By the aid of Eq. (0) and the fact that [fi(X), A*(q)] = 

d>1 Jx ^ ' recurrence relations can be derived 



\»(X),-iP\ 



for the matrix elements of n(X) between the states \<fi n )- 
We have 

{(b m \Xe-^ x \^ n+1 ) = [(n - m - 7) (4>m\ Xe~~> x \</> n ) 
+C m (^ m _!| Xe~^ x \<t> n ) + {4> m \e-i x \<l> n )]/C n+u 
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Then all the necessary matrix elements can be cal- 
culated, starting from terms (4>o\X exp(— jX)\ipo) and 
(00 1 exp(— ~fX) \4>o). These latter can be obtained via 
integration in coordinate representation using the wave 
function </> (x) of Eq. (p|) , and we get 
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(15) 



where if) denotes Euler's digamma function ||. 

The values of the parameters aj, di and ji can be ob- 
tained by fitting the function (|l3|) to experimental dipole 
curves, or to those obtained from molecular calculations. 

The quasi-number states constitute a complete or- 
thonormal basis, hence one can expand the solution 



|$(t)) of Eq. (§ in the form: |*(t)) = V c 



M<t>n)- 



Measuring the time t in units of T = 2tt/ujq, the 
Schrodinger equation can be written as the following in- 
finite system of ordinary differential equations 
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where (0„|i?o|0m) and (4> } are given by the ex- 

pressions ([!(]) and (|l4|), respectively. We note here that 
the usual method of solving the time evolution of the 
molecular state is based on the split operator method 
H . That approach solves the partial differential equa- 
tion corresponding to Eq. (0J , and yields an approximate 
wave function of the variables x and t. In contrast, we 
solve a system of ordinary equations allowing a higher 
accuracy in the calculations. 

We will apply our method to determine the dissocia- 
tion probability of a molecule. This is a time dependent 
quantity defined as the projection probability of |$(t)} 
on Ti + that can be given as 

M M 
V(t) = 1 - £|(<I>(i)hU| 2 = 1 - Y,<n(t)c m (t). (17) 



ro=0 



m=0 



To calculate the dissociation probability above, one has 
to solve the set of the dynamical equations ( |l6| ) and de- 
termine the coefficients c m (t) for < m < [s]. In practice 
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one has to truncate the system (|16j), and we have to use 
a numerical method to find the solution. 

The restriction of the number of the dynamical equa- 
tions to a finite number is done here in two steps. First 
we consider the problem in a large subspace of dimension 
N ^> [s], and neglect the contribution of those states \<j) n ) 
for which n > N . Then the Hamiltonian is represented 
by an operator which is restricted to the finite N di- 
mensional subspace, spanned by the orthonormal system: 
{\if> m ), m = 0,l,2...[s], \<p n ), n=[s] + l,... N}.We 
denote this operator by i?W = Hq — ^ NS) E{i) where 
Hq N ^ and n^E(t) are the truncated operators, so that 
their matrices are of dimension N. 

Before solving the already finite number of ordinary 
differential equations, we perform a unitary transforma- 
tion on our truncated basis and bring Hq into diagonal 
form for n > [s], too. Due to the fact that the matrix of 
Hq is finite tridiagonal this can be done by the help of a 
very fast algorithm. This yields positive energy eigenval- 
ues En (H < n < N) and eigenstates \ipn)i which wc 
can use for the description of the dissociated molecule. 
If N is large enough, then the lowermost positive en- 
ergies En follow densely each other, and approximate 
satisfactorily the continuous energy spectrum above the 
dissociation threshold. In the second step of the approx- 
imation, we restrict the calculation to the bound sub- 
space and to those \il>n) states with n > [s], for which 
the dipole couplings with the bound energy eigenstates 
are non-negligible. In the application to be discussed be- 
low, it turns out that the latter are significant only for 
those states that correspond to the lowest positive energy 
values E^ N) with n = [s] + 1, . . . M, so that M <C N. 
Therefore, in order to follow the dissociation process, it 
is enough to solve Eq. (16) in the basis \i/) m ), 
which means that we expand the time dependent state 
as | <!>(<)) = Xm=o b n {t)\^n), where now the notation \ip n ) 
has been used also for \ip^ ) with n = [s] + 1, . . . M. Then 
we solve the system of M ordinary differential equations 
that describes the time evolution of the coefficients b n 
and calculate the dissociation probability via Eq. jl7| ) 
with the obvious c m (t) — > b m (t) replacement. 

In what follows we are going to apply our method for 
the determination of the dissociation probability of the 
nitrogen-oxide (NO) molecule, excited by an appropri- 
ate laser field H. We have used a Morse potential with 
the following data in the electronic ground state p0[ : 
m = 7.46 a.u., D = 6.497 eV and a = 27.68 nm^, 
yielding s = 54.54. The potential and the dependence 
of the dipole moment on the interatomic separation are 
depicted in Fig. [l| The circles show the dipole moments 
of NO at various interatomic separations that we have 
calculated by an unrestricted density functional method 
(□J. The dotted line is a fitted function of the form 
Eq. ( p"3| ) with two terms, and with the parameter val- 
ues: ai = -9.66, a 2 = 10.64, di = d 2 = 0, 71 = 0.927, 
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FIG. 1: The Morse potential and the dipole moment of the 
NO molecule. 



72 = 0.870. 

The corresponding matrix elements have been deter- 
mined analytically by using the recurrence relations ( [l4| - 
|l5| ), and those between the first second and third neigh- 
bors are shown in Fig ^|. In our calculations we have 
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FIG. 2: Absolute values of the dipole matrix elements be- 
tween first, second, and third neighbor states. The matrix 
elements far from the diagonal are diminishingly small. 

chosen the truncation indices to be N — 3000, and 
M = 200, as increasing these values had no effect on the 
results. (Note that the number of bound states is now 
[s] + 1 = 55). The initial state was taken to be the vibra- 
tional ground state, and we devised the excitation with 
the aim to acheive a significant dissociation probability. 

According to (l2) one can think of the desired time evo- 
lution as a sequence of two-level transitions. There are 
two conditions to be fulfilled: the electric field should be 
approximately resonant with the actual transition, and 
its area per transition 113] should be around 7r. In order 
to maintain the resonance condition during the whole 
process, it is appropriate to use continuously chirped 
pulses Q . If the frequency of the pulse decreases slowly 
enough, then it is possible to consider its effect as a se- 
quence of approximately resonant two- level transitions. 
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For low vibrational quantum numbers the coupling is 
largest between the nearest neighbors, therefore at the 
beginning the excitation has to drive the system through 
the ladder of neighboring energy eigenstates towards the 
continuum. However, as it is seen in Fig. there is a 
minimum in the m — > m + 1 couplings at m — 42, there- 
fore this vibrational quantum number corresponds to a 
trapping state. This is a consequence of the shape of the 
dipole moment curve (Fig. |l|) and it would be absent 
in the linear approximation fJ-(X) ~ X. One can cir- 
cumvent the problem by lengthening the pulse in time, 
but we present a different route here with the three pulses 
shown in Fig. ^| exploiting other than first neighbor tran- 
sitions, as well. The area condition could be met for all 
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FIG. 3: Electric field strength envelopes and time dependent 
frequencies of the applied pulses. Dissociation probability is 
also shown as a function of time measured in units of T — 
2-k/ujq. All three pulses have a distorted secant-hyperbolic 
envelope (see text) with the data: pulse (1): ry = 10.27, 
T t = 246 A, t = 41.0; pulse (2): r f = 71.87, n = 61.61, 
t = 1297.5; pulse (3): r f = 71.87, n = 61.61, to = 2087.6. 
Peak values for all three pulses are: 1.05 x 10 8 V/cm. 

three pulses with distorted secant-hyperbolic envelopes: 
£ = £ a /[exp(io — t)/rf +exp(i — to)/rt] with different ex- 
ponentials in their fronts (jy), and tails (r*). The chirp- 
ing rates are seen in Fig. g, while the other data of the 
pulses are given in the caption. Using these parameters, 
the first pulse drives the molecule into a superposition 
of a few states \tp m ) around m = 31. As seen in Fig. 
||, the m — > m + 2 matrix elements are getting larger 
here than those connecting the neighboring states, and 
the former have also nonvanishing values around m = 42. 
This implies that a second laser pulse with an initial car- 
rier frequency being in approximate resonance with the 
corresponding m — > m+2 transition will continue the dis- 
sociation process and let the molecular state jump over 
the trap. After this second pulse, the distribution of the 
populated states is not so narrow as it was before, which 
is a consequence of the presence of second neighbor tran- 
sitions. The most probably populated states are now 



so close to the continuum limit that a final pulse - even 
without chirping - leads to a dissociation probability that 
is more than 40%. As shown in Fig. ||, this probability 
is practically zero before the third pulse. 

In summary, we have presented a method based on 
supersymmetric quantum mechanics that leads to a sys- 
tem of ordinary differential equations describing molec- 
ular time evolution. We have applied this method for 
devising laser pulses leading to a significant dissociation 
probability of the NO molecule with a realistic dipole 
moment. 
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